library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.2.1 ──
## ✔ ggplot2 3.2.0 ✔ purrr 0.3.2
## ✔ tibble 2.1.1 ✔ dplyr 0.8.3
## ✔ tidyr 1.0.0 ✔ stringr 1.4.0
## ✔ readr 1.3.1 ✔ forcats 0.4.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
library(glue)
##
## Attaching package: 'glue'
## The following object is masked from 'package:dplyr':
##
## collapse
library(tidybayes)
library(cowplot)
##
## ********************************************************
## Note: As of version 1.0.0, cowplot does not change the
## default ggplot2 theme anymore. To recover the previous
## behavior, execute:
## theme_set(theme_cowplot())
## ********************************************************
library(scales)
##
## Attaching package: 'scales'
## The following object is masked from 'package:purrr':
##
## discard
## The following object is masked from 'package:readr':
##
## col_factor
library(GGally)
##
## Attaching package: 'GGally'
## The following object is masked from 'package:dplyr':
##
## nasa
library(knitr)
theme_set(theme_cowplot())
load_existing_sim <- TRUE
Continuing from m1, the most noticeable problem with the initial priors was that too many extremes values were being observed at the subject level, both in terms of cricSD and pMem. This model will shrink the hierarchcial variance parameters for circSD and pMem conditions. It might turn out that the mean priors need to shrunk as well/instead.
\[ \begin{aligned} sigma\_\alpha_0 \sim Normal^+(0, 0.5) ~~ ...~&to~... ~~ sigma\_\alpha_0 \sim Normal^+(0, 0.25) \\ sigma\_\beta_0 \sim Normal^+(0, 1.5) ~~ ...~&to~... ~~ sigma\_\beta_0 \sim Normal^+(0, 0.5) \\ \\ sigma\_\alpha_\Delta \sim Normal^+(0, 0.5) ~~ ...~&to~... ~~ sigma\_\alpha_\Delta \sim Normal^+(0, 0.25) \\ sigma\_\beta_\Delta \sim Normal^+(0, 1) ~~ ...~&to~... ~~ sigma\_\beta_\Delta \sim Normal^+(0, 0.5) \end{aligned} \]
Written down the model is:
\[ \begin{aligned} \mathrm{-likelihood-} \\ error_i &\sim (pMem_i)*VM(0, \kappa_i) + (1 - pMem_i)*Unif(-\pi,\pi) \\ \\ \mathrm{-param~transformation-} \\ \kappa_i &= sd2k(circ\_sd_i) \\ \\ \mathrm{-linear~model-} \\ circ\_sd_i &= exp(\alpha_{0,SUBJ[i]} + \alpha_{\Delta,SUBJ[i]} * postCond) ~~ \mathrm{-log~link~on~sd-} \\ pMem_i &= inv\_logit(\beta_{0,SUBJ[i]} + \beta_{\Delta,SUBJ[i]} * postCond) ~~ \mathrm{-logit~link~on~pMem-} \\ \\ \mathrm{-priors:~all~independent-}\\ \alpha_{0,SUBJ[...]} &\sim Normal(mu\_\alpha_0, sigma\_\alpha_0) \\ \alpha_{\Delta,SUBJ[...]} &\sim Normal(mu\_\alpha_{\Delta}, sigma\_\alpha_{\Delta}) \\ \beta_{0,SUBJ[...]} &\sim Normal(mu\_\beta_0, sigma\_\beta_0) \\ \beta_{\Delta,SUBJ[...]} &\sim Normal(mu\_\beta_{\Delta}, sigma\_\beta_{\Delta}) \\ \\ mu\_\alpha_0 &\sim Normal(4, 0.5)\\ sigma\_\alpha_0 &\sim Normal^+(0, 0.25) \\ mu\_\alpha_{\Delta} &\sim Normal(0, 0.5)\\ sigma\_\alpha_{\Delta} &\sim Normal^+(0, 0.25)\\ mu\_\beta_0 &\sim Normal(0, 1.5)\\ sigma\_\beta_0 &\sim Normal^+(0, 0.5)\\ mu\_\beta_{\Delta} &\sim Normal(0, 1)\\ sigma\_\beta_{\Delta} &\sim Normal^+(0, 0.5)\\ \end{aligned} \]
I am considering the color stimulation data. samples sizes:
obs_data <- read_csv('../data/stimulation_obvs.csv')
## Parsed with column specification:
## cols(
## subj = col_double(),
## subj_index = col_double(),
## stimulation = col_double(),
## error = col_double()
## )
#summarize subj num, obs count, and obs condition split
(subj_summary <- obs_data %>%
group_by(subj_index) %>%
summarise(n_obs = n(), frac_stim = mean(stimulation)))
## # A tibble: 2 x 3
## subj_index n_obs frac_stim
## <dbl> <int> <dbl>
## 1 1 252 0.5
## 2 2 252 0.5
2 subjs with 252 observations each, 126 per condition.
9/12 - I am thinking I should simulate varying effects using the actual samples sizes I have. I also think that this shouldnt matter (at least for a model like this). I’ll also try with a larger number of subjects.
nsubj_sim <- 15
nobs_per_cond_sim <- 126
functions
inv_logit_scalar <- function(x){
return(exp(x)/(exp(x) + 1))
}
inv_logit <- Vectorize(inv_logit_scalar)
sd2k <- function(sd_input) {
#sd_input needs to be radians
#original matlab code
# R <- exp(-S.^2/2);
# K = 1./(R.^3 - 4 * R.^2 + 3 * R);
# K(R < 0.85) = -0.4 + 1.39 * R(R < 0.85) + 0.43./(1 - R(R < 0.85));
# K(R < 0.53) = 2 * R(R < 0.53) + R(R < 0.53).^3 + (5 * R(R < 0.53).^5)/6;
R <- exp(-(sd_input^2)/2)
K <- 1/((R^3) - (4 * R^2) + (3*R))
if(R < 0.85){
K <- -0.4 + (1.39 * R) + (0.43/(1-R))
}
if(R < 0.53){
K <- (2 * R) + R^3 + (5 * R^5)/6
}
return(K)
}
sd2k_vec <- Vectorize(sd2k)
# this function takes a single set of parameters and covariates and simlates an observation from the mixture model
simulateData_likelihood <- function(pMem,
k,
nobs = 1){
# which part of the mixture does this obs come from
memFlip <- rbernoulli(nobs, pMem)
obs <- memFlip * ( Rfast::rvonmises(nobs, pi, k) - pi) + (1 - memFlip) * runif(nobs, -pi, pi)
obs <- as.numeric(obs)
obs_tibble <- tibble(obs_radian = obs, obs_degree = pracma::rad2deg(obs))
return(obs_tibble)
}
simulateData_subj <- function(subj_alpha0, subj_alphaD,
subj_beta0, subj_betaD,
nobs_per_condition){
conditions <- c(0,1)
stimulation <- rep(conditions, each = nobs_per_condition)
params <- list(pMem = inv_logit(subj_beta0 + (subj_betaD * stimulation)),
k = sd2k_vec(
pracma::deg2rad(
exp(subj_alpha0 + (subj_alphaD * stimulation)))))
subj_obs <- pmap_dfr(params, simulateData_likelihood) %>%
cbind(stimulation)
return(subj_obs)
}
# this function will take a single set of simulated parameter draws and then simulate an implied data set
simulateData <- function(sim_num,
alpha0_mu, alpha0_sigma,
alphaD_mu, alphaD_sigma,
beta0_mu, beta0_sigma,
betaD_mu, betaD_sigma,
nsubj = 100, nobs_per_condition = 100){
print(glue("simulation: {sim_num}"))
# generate lower level subject-specific parameters
sim_subj_alpha0 <- rnorm(nsubj, alpha0_mu, alpha0_sigma)
sim_subj_alphaD <- rnorm(nsubj, alphaD_mu, alphaD_sigma)
sim_subj_beta0 <- rnorm(nsubj, beta0_mu, beta0_sigma)
sim_subj_betaD <- rnorm(nsubj, betaD_mu, betaD_sigma)
sim_data <- tibble(
subj = 1:nsubj,
subj_alpha0 = sim_subj_alpha0,
subj_alphaD = sim_subj_alphaD,
subj_beta0 = sim_subj_beta0,
subj_betaD = sim_subj_betaD,
nobs_per_condition = nobs_per_condition
)
sim_data <- sim_data %>%
mutate(subj_obs = pmap(sim_data %>% select(-subj), simulateData_subj))
return(sim_data)
}
run simulate
found_existing_sim <- FALSE
if (file.exists("sim_datasets.rds")){
found_existing_sim <- TRUE
}
if (load_existing_sim == FALSE || found_existing_sim == FALSE){
nsim_datasets <- 1200
sim_priors <- tibble(
sim_num = 1:nsim_datasets,
alpha0_mu = rnorm(nsim_datasets, 4, 0.5),
alpha0_sigma = abs(rnorm(nsim_datasets, 0, 0.25)),
alphaD_mu = rnorm(nsim_datasets, 0, 0.5),
alphaD_sigma = abs(rnorm(nsim_datasets, 0, 0.25)),
beta0_mu = rnorm(nsim_datasets, 0, 1.5),
beta0_sigma = abs(rnorm(nsim_datasets, 0, 0.5)),
betaD_mu = rnorm(nsim_datasets, 0, 1),
betaD_sigma = abs(rnorm(nsim_datasets, 0, 0.5)),
nsubj = nsubj_sim,
nobs_per_cond = nobs_per_cond_sim
)
sim_datasets <- sim_priors %>%
mutate(dataset = pmap(sim_priors, simulateData))
saveRDS(sim_datasets, file = "sim_datasets.rds")
}else{
sim_datasets <- readRDS("sim_datasets.rds")
}
## simulation: 1
## simulation: 2
## simulation: 3
## simulation: 4
## simulation: 5
## simulation: 6
## simulation: 7
## simulation: 8
## simulation: 9
## simulation: 10
## simulation: 11
## simulation: 12
## simulation: 13
## simulation: 14
## simulation: 15
## simulation: 16
## simulation: 17
## simulation: 18
## simulation: 19
## simulation: 20
## simulation: 21
## simulation: 22
## simulation: 23
## simulation: 24
## simulation: 25
## simulation: 26
## simulation: 27
## simulation: 28
## simulation: 29
## simulation: 30
## simulation: 31
## simulation: 32
## simulation: 33
## simulation: 34
## simulation: 35
## simulation: 36
## simulation: 37
## simulation: 38
## simulation: 39
## simulation: 40
## simulation: 41
## simulation: 42
## simulation: 43
## simulation: 44
## simulation: 45
## simulation: 46
## simulation: 47
## simulation: 48
## simulation: 49
## simulation: 50
## simulation: 51
## simulation: 52
## simulation: 53
## simulation: 54
## simulation: 55
## simulation: 56
## simulation: 57
## simulation: 58
## simulation: 59
## simulation: 60
## simulation: 61
## simulation: 62
## simulation: 63
## simulation: 64
## simulation: 65
## simulation: 66
## simulation: 67
## simulation: 68
## simulation: 69
## simulation: 70
## simulation: 71
## simulation: 72
## simulation: 73
## simulation: 74
## simulation: 75
## simulation: 76
## simulation: 77
## simulation: 78
## simulation: 79
## simulation: 80
## simulation: 81
## simulation: 82
## simulation: 83
## simulation: 84
## simulation: 85
## simulation: 86
## simulation: 87
## simulation: 88
## simulation: 89
## simulation: 90
## simulation: 91
## simulation: 92
## simulation: 93
## simulation: 94
## simulation: 95
## simulation: 96
## simulation: 97
## simulation: 98
## simulation: 99
## simulation: 100
## simulation: 101
## simulation: 102
## simulation: 103
## simulation: 104
## simulation: 105
## simulation: 106
## simulation: 107
## simulation: 108
## simulation: 109
## simulation: 110
## simulation: 111
## simulation: 112
## simulation: 113
## simulation: 114
## simulation: 115
## simulation: 116
## simulation: 117
## simulation: 118
## simulation: 119
## simulation: 120
## simulation: 121
## simulation: 122
## simulation: 123
## simulation: 124
## simulation: 125
## simulation: 126
## simulation: 127
## simulation: 128
## simulation: 129
## simulation: 130
## simulation: 131
## simulation: 132
## simulation: 133
## simulation: 134
## simulation: 135
## simulation: 136
## simulation: 137
## simulation: 138
## simulation: 139
## simulation: 140
## simulation: 141
## simulation: 142
## simulation: 143
## simulation: 144
## simulation: 145
## simulation: 146
## simulation: 147
## simulation: 148
## simulation: 149
## simulation: 150
## simulation: 151
## simulation: 152
## simulation: 153
## simulation: 154
## simulation: 155
## simulation: 156
## simulation: 157
## simulation: 158
## simulation: 159
## simulation: 160
## simulation: 161
## simulation: 162
## simulation: 163
## simulation: 164
## simulation: 165
## simulation: 166
## simulation: 167
## simulation: 168
## simulation: 169
## simulation: 170
## simulation: 171
## simulation: 172
## simulation: 173
## simulation: 174
## simulation: 175
## simulation: 176
## simulation: 177
## simulation: 178
## simulation: 179
## simulation: 180
## simulation: 181
## simulation: 182
## simulation: 183
## simulation: 184
## simulation: 185
## simulation: 186
## simulation: 187
## simulation: 188
## simulation: 189
## simulation: 190
## simulation: 191
## simulation: 192
## simulation: 193
## simulation: 194
## simulation: 195
## simulation: 196
## simulation: 197
## simulation: 198
## simulation: 199
## simulation: 200
## simulation: 201
## simulation: 202
## simulation: 203
## simulation: 204
## simulation: 205
## simulation: 206
## simulation: 207
## simulation: 208
## simulation: 209
## simulation: 210
## simulation: 211
## simulation: 212
## simulation: 213
## simulation: 214
## simulation: 215
## simulation: 216
## simulation: 217
## simulation: 218
## simulation: 219
## simulation: 220
## simulation: 221
## simulation: 222
## simulation: 223
## simulation: 224
## simulation: 225
## simulation: 226
## simulation: 227
## simulation: 228
## simulation: 229
## simulation: 230
## simulation: 231
## simulation: 232
## simulation: 233
## simulation: 234
## simulation: 235
## simulation: 236
## simulation: 237
## simulation: 238
## simulation: 239
## simulation: 240
## simulation: 241
## simulation: 242
## simulation: 243
## simulation: 244
## simulation: 245
## simulation: 246
## simulation: 247
## simulation: 248
## simulation: 249
## simulation: 250
## simulation: 251
## simulation: 252
## simulation: 253
## simulation: 254
## simulation: 255
## simulation: 256
## simulation: 257
## simulation: 258
## simulation: 259
## simulation: 260
## simulation: 261
## simulation: 262
## simulation: 263
## simulation: 264
## simulation: 265
## simulation: 266
## simulation: 267
## simulation: 268
## simulation: 269
## simulation: 270
## simulation: 271
## simulation: 272
## simulation: 273
## simulation: 274
## simulation: 275
## simulation: 276
## simulation: 277
## simulation: 278
## simulation: 279
## simulation: 280
## simulation: 281
## simulation: 282
## simulation: 283
## simulation: 284
## simulation: 285
## simulation: 286
## simulation: 287
## simulation: 288
## simulation: 289
## simulation: 290
## simulation: 291
## simulation: 292
## simulation: 293
## simulation: 294
## simulation: 295
## simulation: 296
## simulation: 297
## simulation: 298
## simulation: 299
## simulation: 300
## simulation: 301
## simulation: 302
## simulation: 303
## simulation: 304
## simulation: 305
## simulation: 306
## simulation: 307
## simulation: 308
## simulation: 309
## simulation: 310
## simulation: 311
## simulation: 312
## simulation: 313
## simulation: 314
## simulation: 315
## simulation: 316
## simulation: 317
## simulation: 318
## simulation: 319
## simulation: 320
## simulation: 321
## simulation: 322
## simulation: 323
## simulation: 324
## simulation: 325
## simulation: 326
## simulation: 327
## simulation: 328
## simulation: 329
## simulation: 330
## simulation: 331
## simulation: 332
## simulation: 333
## simulation: 334
## simulation: 335
## simulation: 336
## simulation: 337
## simulation: 338
## simulation: 339
## simulation: 340
## simulation: 341
## simulation: 342
## simulation: 343
## simulation: 344
## simulation: 345
## simulation: 346
## simulation: 347
## simulation: 348
## simulation: 349
## simulation: 350
## simulation: 351
## simulation: 352
## simulation: 353
## simulation: 354
## simulation: 355
## simulation: 356
## simulation: 357
## simulation: 358
## simulation: 359
## simulation: 360
## simulation: 361
## simulation: 362
## simulation: 363
## simulation: 364
## simulation: 365
## simulation: 366
## simulation: 367
## simulation: 368
## simulation: 369
## simulation: 370
## simulation: 371
## simulation: 372
## simulation: 373
## simulation: 374
## simulation: 375
## simulation: 376
## simulation: 377
## simulation: 378
## simulation: 379
## simulation: 380
## simulation: 381
## simulation: 382
## simulation: 383
## simulation: 384
## simulation: 385
## simulation: 386
## simulation: 387
## simulation: 388
## simulation: 389
## simulation: 390
## simulation: 391
## simulation: 392
## simulation: 393
## simulation: 394
## simulation: 395
## simulation: 396
## simulation: 397
## simulation: 398
## simulation: 399
## simulation: 400
## simulation: 401
## simulation: 402
## simulation: 403
## simulation: 404
## simulation: 405
## simulation: 406
## simulation: 407
## simulation: 408
## simulation: 409
## simulation: 410
## simulation: 411
## simulation: 412
## simulation: 413
## simulation: 414
## simulation: 415
## simulation: 416
## simulation: 417
## simulation: 418
## simulation: 419
## simulation: 420
## simulation: 421
## simulation: 422
## simulation: 423
## simulation: 424
## simulation: 425
## simulation: 426
## simulation: 427
## simulation: 428
## simulation: 429
## simulation: 430
## simulation: 431
## simulation: 432
## simulation: 433
## simulation: 434
## simulation: 435
## simulation: 436
## simulation: 437
## simulation: 438
## simulation: 439
## simulation: 440
## simulation: 441
## simulation: 442
## simulation: 443
## simulation: 444
## simulation: 445
## simulation: 446
## simulation: 447
## simulation: 448
## simulation: 449
## simulation: 450
## simulation: 451
## simulation: 452
## simulation: 453
## simulation: 454
## simulation: 455
## simulation: 456
## simulation: 457
## simulation: 458
## simulation: 459
## simulation: 460
## simulation: 461
## simulation: 462
## simulation: 463
## simulation: 464
## simulation: 465
## simulation: 466
## simulation: 467
## simulation: 468
## simulation: 469
## simulation: 470
## simulation: 471
## simulation: 472
## simulation: 473
## simulation: 474
## simulation: 475
## simulation: 476
## simulation: 477
## simulation: 478
## simulation: 479
## simulation: 480
## simulation: 481
## simulation: 482
## simulation: 483
## simulation: 484
## simulation: 485
## simulation: 486
## simulation: 487
## simulation: 488
## simulation: 489
## simulation: 490
## simulation: 491
## simulation: 492
## simulation: 493
## simulation: 494
## simulation: 495
## simulation: 496
## simulation: 497
## simulation: 498
## simulation: 499
## simulation: 500
## simulation: 501
## simulation: 502
## simulation: 503
## simulation: 504
## simulation: 505
## simulation: 506
## simulation: 507
## simulation: 508
## simulation: 509
## simulation: 510
## simulation: 511
## simulation: 512
## simulation: 513
## simulation: 514
## simulation: 515
## simulation: 516
## simulation: 517
## simulation: 518
## simulation: 519
## simulation: 520
## simulation: 521
## simulation: 522
## simulation: 523
## simulation: 524
## simulation: 525
## simulation: 526
## simulation: 527
## simulation: 528
## simulation: 529
## simulation: 530
## simulation: 531
## simulation: 532
## simulation: 533
## simulation: 534
## simulation: 535
## simulation: 536
## simulation: 537
## simulation: 538
## simulation: 539
## simulation: 540
## simulation: 541
## simulation: 542
## simulation: 543
## simulation: 544
## simulation: 545
## simulation: 546
## simulation: 547
## simulation: 548
## simulation: 549
## simulation: 550
## simulation: 551
## simulation: 552
## simulation: 553
## simulation: 554
## simulation: 555
## simulation: 556
## simulation: 557
## simulation: 558
## simulation: 559
## simulation: 560
## simulation: 561
## simulation: 562
## simulation: 563
## simulation: 564
## simulation: 565
## simulation: 566
## simulation: 567
## simulation: 568
## simulation: 569
## simulation: 570
## simulation: 571
## simulation: 572
## simulation: 573
## simulation: 574
## simulation: 575
## simulation: 576
## simulation: 577
## simulation: 578
## simulation: 579
## simulation: 580
## simulation: 581
## simulation: 582
## simulation: 583
## simulation: 584
## simulation: 585
## simulation: 586
## simulation: 587
## simulation: 588
## simulation: 589
## simulation: 590
## simulation: 591
## simulation: 592
## simulation: 593
## simulation: 594
## simulation: 595
## simulation: 596
## simulation: 597
## simulation: 598
## simulation: 599
## simulation: 600
## simulation: 601
## simulation: 602
## simulation: 603
## simulation: 604
## simulation: 605
## simulation: 606
## simulation: 607
## simulation: 608
## simulation: 609
## simulation: 610
## simulation: 611
## simulation: 612
## simulation: 613
## simulation: 614
## simulation: 615
## simulation: 616
## simulation: 617
## simulation: 618
## simulation: 619
## simulation: 620
## simulation: 621
## simulation: 622
## simulation: 623
## simulation: 624
## simulation: 625
## simulation: 626
## simulation: 627
## simulation: 628
## simulation: 629
## simulation: 630
## simulation: 631
## simulation: 632
## simulation: 633
## simulation: 634
## simulation: 635
## simulation: 636
## simulation: 637
## simulation: 638
## simulation: 639
## simulation: 640
## simulation: 641
## simulation: 642
## simulation: 643
## simulation: 644
## simulation: 645
## simulation: 646
## simulation: 647
## simulation: 648
## simulation: 649
## simulation: 650
## simulation: 651
## simulation: 652
## simulation: 653
## simulation: 654
## simulation: 655
## simulation: 656
## simulation: 657
## simulation: 658
## simulation: 659
## simulation: 660
## simulation: 661
## simulation: 662
## simulation: 663
## simulation: 664
## simulation: 665
## simulation: 666
## simulation: 667
## simulation: 668
## simulation: 669
## simulation: 670
## simulation: 671
## simulation: 672
## simulation: 673
## simulation: 674
## simulation: 675
## simulation: 676
## simulation: 677
## simulation: 678
## simulation: 679
## simulation: 680
## simulation: 681
## simulation: 682
## simulation: 683
## simulation: 684
## simulation: 685
## simulation: 686
## simulation: 687
## simulation: 688
## simulation: 689
## simulation: 690
## simulation: 691
## simulation: 692
## simulation: 693
## simulation: 694
## simulation: 695
## simulation: 696
## simulation: 697
## simulation: 698
## simulation: 699
## simulation: 700
## simulation: 701
## simulation: 702
## simulation: 703
## simulation: 704
## simulation: 705
## simulation: 706
## simulation: 707
## simulation: 708
## simulation: 709
## simulation: 710
## simulation: 711
## simulation: 712
## simulation: 713
## simulation: 714
## simulation: 715
## simulation: 716
## simulation: 717
## simulation: 718
## simulation: 719
## simulation: 720
## simulation: 721
## simulation: 722
## simulation: 723
## simulation: 724
## simulation: 725
## simulation: 726
## simulation: 727
## simulation: 728
## simulation: 729
## simulation: 730
## simulation: 731
## simulation: 732
## simulation: 733
## simulation: 734
## simulation: 735
## simulation: 736
## simulation: 737
## simulation: 738
## simulation: 739
## simulation: 740
## simulation: 741
## simulation: 742
## simulation: 743
## simulation: 744
## simulation: 745
## simulation: 746
## simulation: 747
## simulation: 748
## simulation: 749
## simulation: 750
## simulation: 751
## simulation: 752
## simulation: 753
## simulation: 754
## simulation: 755
## simulation: 756
## simulation: 757
## simulation: 758
## simulation: 759
## simulation: 760
## simulation: 761
## simulation: 762
## simulation: 763
## simulation: 764
## simulation: 765
## simulation: 766
## simulation: 767
## simulation: 768
## simulation: 769
## simulation: 770
## simulation: 771
## simulation: 772
## simulation: 773
## simulation: 774
## simulation: 775
## simulation: 776
## simulation: 777
## simulation: 778
## simulation: 779
## simulation: 780
## simulation: 781
## simulation: 782
## simulation: 783
## simulation: 784
## simulation: 785
## simulation: 786
## simulation: 787
## simulation: 788
## simulation: 789
## simulation: 790
## simulation: 791
## simulation: 792
## simulation: 793
## simulation: 794
## simulation: 795
## simulation: 796
## simulation: 797
## simulation: 798
## simulation: 799
## simulation: 800
## simulation: 801
## simulation: 802
## simulation: 803
## simulation: 804
## simulation: 805
## simulation: 806
## simulation: 807
## simulation: 808
## simulation: 809
## simulation: 810
## simulation: 811
## simulation: 812
## simulation: 813
## simulation: 814
## simulation: 815
## simulation: 816
## simulation: 817
## simulation: 818
## simulation: 819
## simulation: 820
## simulation: 821
## simulation: 822
## simulation: 823
## simulation: 824
## simulation: 825
## simulation: 826
## simulation: 827
## simulation: 828
## simulation: 829
## simulation: 830
## simulation: 831
## simulation: 832
## simulation: 833
## simulation: 834
## simulation: 835
## simulation: 836
## simulation: 837
## simulation: 838
## simulation: 839
## simulation: 840
## simulation: 841
## simulation: 842
## simulation: 843
## simulation: 844
## simulation: 845
## simulation: 846
## simulation: 847
## simulation: 848
## simulation: 849
## simulation: 850
## simulation: 851
## simulation: 852
## simulation: 853
## simulation: 854
## simulation: 855
## simulation: 856
## simulation: 857
## simulation: 858
## simulation: 859
## simulation: 860
## simulation: 861
## simulation: 862
## simulation: 863
## simulation: 864
## simulation: 865
## simulation: 866
## simulation: 867
## simulation: 868
## simulation: 869
## simulation: 870
## simulation: 871
## simulation: 872
## simulation: 873
## simulation: 874
## simulation: 875
## simulation: 876
## simulation: 877
## simulation: 878
## simulation: 879
## simulation: 880
## simulation: 881
## simulation: 882
## simulation: 883
## simulation: 884
## simulation: 885
## simulation: 886
## simulation: 887
## simulation: 888
## simulation: 889
## simulation: 890
## simulation: 891
## simulation: 892
## simulation: 893
## simulation: 894
## simulation: 895
## simulation: 896
## simulation: 897
## simulation: 898
## simulation: 899
## simulation: 900
## simulation: 901
## simulation: 902
## simulation: 903
## simulation: 904
## simulation: 905
## simulation: 906
## simulation: 907
## simulation: 908
## simulation: 909
## simulation: 910
## simulation: 911
## simulation: 912
## simulation: 913
## simulation: 914
## simulation: 915
## simulation: 916
## simulation: 917
## simulation: 918
## simulation: 919
## simulation: 920
## simulation: 921
## simulation: 922
## simulation: 923
## simulation: 924
## simulation: 925
## simulation: 926
## simulation: 927
## simulation: 928
## simulation: 929
## simulation: 930
## simulation: 931
## simulation: 932
## simulation: 933
## simulation: 934
## simulation: 935
## simulation: 936
## simulation: 937
## simulation: 938
## simulation: 939
## simulation: 940
## simulation: 941
## simulation: 942
## simulation: 943
## simulation: 944
## simulation: 945
## simulation: 946
## simulation: 947
## simulation: 948
## simulation: 949
## simulation: 950
## simulation: 951
## simulation: 952
## simulation: 953
## simulation: 954
## simulation: 955
## simulation: 956
## simulation: 957
## simulation: 958
## simulation: 959
## simulation: 960
## simulation: 961
## simulation: 962
## simulation: 963
## simulation: 964
## simulation: 965
## simulation: 966
## simulation: 967
## simulation: 968
## simulation: 969
## simulation: 970
## simulation: 971
## simulation: 972
## simulation: 973
## simulation: 974
## simulation: 975
## simulation: 976
## simulation: 977
## simulation: 978
## simulation: 979
## simulation: 980
## simulation: 981
## simulation: 982
## simulation: 983
## simulation: 984
## simulation: 985
## simulation: 986
## simulation: 987
## simulation: 988
## simulation: 989
## simulation: 990
## simulation: 991
## simulation: 992
## simulation: 993
## simulation: 994
## simulation: 995
## simulation: 996
## simulation: 997
## simulation: 998
## simulation: 999
## simulation: 1000
## simulation: 1001
## simulation: 1002
## simulation: 1003
## simulation: 1004
## simulation: 1005
## simulation: 1006
## simulation: 1007
## simulation: 1008
## simulation: 1009
## simulation: 1010
## simulation: 1011
## simulation: 1012
## simulation: 1013
## simulation: 1014
## simulation: 1015
## simulation: 1016
## simulation: 1017
## simulation: 1018
## simulation: 1019
## simulation: 1020
## simulation: 1021
## simulation: 1022
## simulation: 1023
## simulation: 1024
## simulation: 1025
## simulation: 1026
## simulation: 1027
## simulation: 1028
## simulation: 1029
## simulation: 1030
## simulation: 1031
## simulation: 1032
## simulation: 1033
## simulation: 1034
## simulation: 1035
## simulation: 1036
## simulation: 1037
## simulation: 1038
## simulation: 1039
## simulation: 1040
## simulation: 1041
## simulation: 1042
## simulation: 1043
## simulation: 1044
## simulation: 1045
## simulation: 1046
## simulation: 1047
## simulation: 1048
## simulation: 1049
## simulation: 1050
## simulation: 1051
## simulation: 1052
## simulation: 1053
## simulation: 1054
## simulation: 1055
## simulation: 1056
## simulation: 1057
## simulation: 1058
## simulation: 1059
## simulation: 1060
## simulation: 1061
## simulation: 1062
## simulation: 1063
## simulation: 1064
## simulation: 1065
## simulation: 1066
## simulation: 1067
## simulation: 1068
## simulation: 1069
## simulation: 1070
## simulation: 1071
## simulation: 1072
## simulation: 1073
## simulation: 1074
## simulation: 1075
## simulation: 1076
## simulation: 1077
## simulation: 1078
## simulation: 1079
## simulation: 1080
## simulation: 1081
## simulation: 1082
## simulation: 1083
## simulation: 1084
## simulation: 1085
## simulation: 1086
## simulation: 1087
## simulation: 1088
## simulation: 1089
## simulation: 1090
## simulation: 1091
## simulation: 1092
## simulation: 1093
## simulation: 1094
## simulation: 1095
## simulation: 1096
## simulation: 1097
## simulation: 1098
## simulation: 1099
## simulation: 1100
## simulation: 1101
## simulation: 1102
## simulation: 1103
## simulation: 1104
## simulation: 1105
## simulation: 1106
## simulation: 1107
## simulation: 1108
## simulation: 1109
## simulation: 1110
## simulation: 1111
## simulation: 1112
## simulation: 1113
## simulation: 1114
## simulation: 1115
## simulation: 1116
## simulation: 1117
## simulation: 1118
## simulation: 1119
## simulation: 1120
## simulation: 1121
## simulation: 1122
## simulation: 1123
## simulation: 1124
## simulation: 1125
## simulation: 1126
## simulation: 1127
## simulation: 1128
## simulation: 1129
## simulation: 1130
## simulation: 1131
## simulation: 1132
## simulation: 1133
## simulation: 1134
## simulation: 1135
## simulation: 1136
## simulation: 1137
## simulation: 1138
## simulation: 1139
## simulation: 1140
## simulation: 1141
## simulation: 1142
## simulation: 1143
## simulation: 1144
## simulation: 1145
## simulation: 1146
## simulation: 1147
## simulation: 1148
## simulation: 1149
## simulation: 1150
## simulation: 1151
## simulation: 1152
## simulation: 1153
## simulation: 1154
## simulation: 1155
## simulation: 1156
## simulation: 1157
## simulation: 1158
## simulation: 1159
## simulation: 1160
## simulation: 1161
## simulation: 1162
## simulation: 1163
## simulation: 1164
## simulation: 1165
## simulation: 1166
## simulation: 1167
## simulation: 1168
## simulation: 1169
## simulation: 1170
## simulation: 1171
## simulation: 1172
## simulation: 1173
## simulation: 1174
## simulation: 1175
## simulation: 1176
## simulation: 1177
## simulation: 1178
## simulation: 1179
## simulation: 1180
## simulation: 1181
## simulation: 1182
## simulation: 1183
## simulation: 1184
## simulation: 1185
## simulation: 1186
## simulation: 1187
## simulation: 1188
## simulation: 1189
## simulation: 1190
## simulation: 1191
## simulation: 1192
## simulation: 1193
## simulation: 1194
## simulation: 1195
## simulation: 1196
## simulation: 1197
## simulation: 1198
## simulation: 1199
## simulation: 1200
Check sim data
head(sim_datasets)
## # A tibble: 6 x 12
## sim_num alpha0_mu alpha0_sigma alphaD_mu alphaD_sigma beta0_mu
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2.74 0.0679 -0.823 0.0797 -2.04
## 2 2 3.63 0.123 1.70 0.414 1.52
## 3 3 3.53 0.214 -0.0641 0.314 0.479
## 4 4 4.32 0.169 1.55 0.0785 -0.110
## 5 5 3.53 0.294 0.183 0.170 0.324
## 6 6 4.41 0.424 0.610 0.358 -0.00737
## # … with 6 more variables: beta0_sigma <dbl>, betaD_mu <dbl>,
## # betaD_sigma <dbl>, nsubj <dbl>, nobs_per_cond <dbl>, dataset <list>
head(sim_datasets$dataset[[1]])
## # A tibble: 6 x 7
## subj subj_alpha0 subj_alphaD subj_beta0 subj_betaD nobs_per_condit…
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 1 2.73 -0.893 -2.04 -0.491 126
## 2 2 2.75 -0.831 -1.91 -0.309 126
## 3 3 2.76 -0.925 -1.98 -0.323 126
## 4 4 2.65 -0.865 -1.91 -0.200 126
## 5 5 2.75 -0.862 -1.67 -0.599 126
## 6 6 2.73 -0.734 -1.73 -0.625 126
## # … with 1 more variable: subj_obs <list>
head(sim_datasets$dataset[[1]]$subj_obs[[1]])
## obs_radian obs_degree stimulation
## 1 2.04052586 116.913520 0
## 2 0.80937892 46.373996 0
## 3 0.08309119 4.760774 0
## 4 -1.93824370 -111.053184 0
## 5 -1.60348221 -91.872763 0
## 6 1.22558395 70.220788 0
Ideas:
Definitely plot the priors on the intercept and slope means, transformed back into sd
plot the prior predicitive densities for subjects
Some other plot ideas:
shaded histogram plot, combining data across stimulation conditions - check if data looks too peaked or too flat
shaded histogram plot, one for each stimulation condition - same check ^
plot of density lines, one from each sim dataset, one with both conditions and one with them split out - this would help identify weird distribution more that shaded histogram perhaps
sim_datasets %>%
select(alpha0_mu, alpha0_sigma, alphaD_mu, alphaD_sigma) %>%
ggpairs(progress = FALSE)
Good.
sim_datasets %>%
select(beta0_mu, beta0_sigma, betaD_mu, betaD_sigma) %>%
ggpairs(progress = FALSE)
Good.
params <- cross_df(list(pMem = seq(0, 1, 0.1), sd = seq(5, 18, 3)))
params %>%
mutate(sims = map2(params$pMem, sd2k_vec(pracma::deg2rad(params$sd)), ~simulateData_likelihood(.x, .y, 1000))) %>%
unnest(sims) %>%
ggplot(aes(x = obs_degree)) +
geom_histogram(binwidth = 1) +
facet_grid(rows = vars(pMem), cols = vars(sd), labeller = 'label_both', scales = 'free')
params <- cross_df(list(pMem = seq(0, 1, 0.1), sd = seq(20, 50, 5)))
params %>%
mutate(sims = map2(params$pMem, sd2k_vec(pracma::deg2rad(params$sd)), ~simulateData_likelihood(.x, .y, 1000))) %>%
unnest(sims) %>%
ggplot(aes(x = obs_degree)) +
geom_histogram(binwidth = 1) +
facet_grid(rows = vars(pMem), cols = vars(sd), labeller = 'label_both', scales = 'free')
params <- cross_df(list(pMem = seq(0, 1, 0.1), sd = seq(55, 105, 5)))
params %>%
mutate(sims = map2(params$pMem, sd2k_vec(pracma::deg2rad(params$sd)), ~simulateData_likelihood(.x, .y, 1000))) %>%
unnest(sims) %>%
ggplot(aes(x = obs_degree)) +
geom_histogram(binwidth = 1)+
facet_grid(rows = vars(pMem), cols = vars(sd), labeller = 'label_both')
(pre group mean, post group mean, post-pre mean change)
No accounting for subject variance here
##
## alpha0_mu
##
## linked prior distribution
alpha0_mu_hist <- sim_datasets %>%
ggplot(aes(x = alpha0_mu)) +
geom_histogram(binwidth = 0.05, fill = "#F16C66") +
theme(legend.position = "none",
axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
labs(x = "linked prior on group mean for circSD, pre condition",
title = "circSD group mean prior, alpha0_mu") +
scale_x_continuous(limits = c(1,7))
##
## alpha0_mu + alphaD_mu = alpha1_mu
##
## linked prior distribution
alpha1_mu_hist <- sim_datasets %>%
transmute(alpha1_mu = alpha0_mu + alphaD_mu) %>%
ggplot(aes(x = alpha1_mu)) +
geom_histogram(binwidth = 0.05, fill = "#685369") +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
labs(x = "linked prior on group mean for circSD, post condition",
title = "circSD group mean prior, alpha0_mu + alphaD_mu") +
scale_x_continuous(limits = c(1,7))
plot_grid(alpha0_mu_hist, alpha1_mu_hist, nrow = 1, align = "h")
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
##
## alpha0_mu
##
## delinked prior alpha0_mu
alpha0_mu_delinked_stats <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu)) %>%
summarise(less_than_5 = sum(delinked < 5)/n(), greater_than_120 = sum(delinked > 120)/n())
alpha0_mu_delinked_hist <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 2, fill = "#F16C66") +
geom_vline(xintercept = c(5, 120), linetype = "dashed") +
labs(x = "prior on group mean for circSD, pre condition",
title = "circSD group mean prior, exp(alpha0_mu)",
subtitle = glue("prob(circSD < 5) = {alpha0_mu_delinked_stats$less_than_5}\nprob(circSD > 120) = {alpha0_mu_delinked_stats$greater_than_120}")) +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12))
##
## alpha0_mu + alphaD_mu = alpha1_mu
##
## delinked prior alpha0_mu + alphaD_mu
alpha1_mu_delinked_hist <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu + alphaD_mu)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 2, fill = "#685369") +
labs(x = "prior on group mean for circSD, post condition",
title = "circSD group mean prior, exp(alpha0_mu + alphaD_mu)") +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12))
plot_grid(alpha0_mu_delinked_hist, alpha1_mu_delinked_hist, nrow = 1, align = "h")
Should actually switch upper bound to 100.
##
## exp(alpha0_mu + alphaD_mu) - exp(alpha0_mu) plot
##
alphaD_mu_delinked_stats <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu + alphaD_mu) - exp(alpha0_mu)) %>%
summarise(less_than_n100 = sum(delinked < -100)/n(), greater_than_100 = sum(delinked > 100)/n())
alphaD_mu_delinked_hist <- sim_datasets %>%
transmute(delinked = exp(alpha0_mu + alphaD_mu) - exp(alpha0_mu)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 2, fill = "#00AEB2") +
geom_vline(xintercept = c(-100, 100), linetype = "dashed") +
labs(x = "prior on group mean for delta-circSD, mean post minus mean pre",
title = "delta-circSD group mean prior, exp(alpha0_mu + alphaD_mu) - exp(alpha0_mu)",
subtitle = glue("prob(delta-circSD < -100) = {alphaD_mu_delinked_stats$less_than_n100}\nprob(delta-circSD > 100) = {alphaD_mu_delinked_stats$greater_than_100}")) +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
scale_x_continuous(breaks= pretty_breaks(10))
alphaD_mu_delinked_hist
Same as m1, still good on mean.
(based on simulated subjects from each dataset)
These plots indicate the effects of the priors on alpha0_mu, alpha0_sigma, alphaD_mu, alphaD_sigma.
(ignoring simulation sample sizes)
sim_datasets_unnest <- sim_datasets %>%
unnest(dataset) %>%
mutate(alpha0_mu_delinked = exp(alpha0_mu),
alpha1_mu_delinked = exp(alpha0_mu + alphaD_mu),
subj_pre_circSD = exp(subj_alpha0),
subj_post_circSD = exp(subj_alphaD + subj_alpha0),
subj_circSD_ES = subj_post_circSD - subj_pre_circSD)
subj_pre_circSD_plot <- sim_datasets_unnest %>%
mutate(subj_pre_circSD = if_else(subj_pre_circSD > 120, 120, subj_pre_circSD)) %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_pre_circSD)) +
geom_histogram(binwidth = 5, fill = "red", alpha = 0.3) +
labs(x = "subj_pre_circSD [prior predictive distribution]",
subtitle = glue("p(<5) = {sum(sim_datasets_unnest$subj_pre_circSD < 5)/length(sim_datasets_unnest$subj_pre_circSD)}\n p(>120) = {sum(sim_datasets_unnest$subj_pre_circSD > 120)/length(sim_datasets_unnest$subj_pre_circSD)}"))
subj_post_circSD_plot <- sim_datasets_unnest %>%
mutate(subj_post_circSD = if_else(subj_post_circSD > 120, 120, subj_post_circSD)) %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_post_circSD)) +
geom_histogram(binwidth = 5, fill = "red", alpha = 0.3) +
labs(x = "subj_post_circSD [prior predictive distribution]",
subtitle = glue("p(<5) = {sum(sim_datasets_unnest$subj_post_circSD < 5)/length(sim_datasets_unnest$subj_post_circSD)}\n p(>120) = {sum(sim_datasets_unnest$subj_post_circSD > 120)/length(sim_datasets_unnest$subj_post_circSD)}"))
subj_circSD_ES_plot <- sim_datasets_unnest %>%
mutate(subj_circSD_ES = if_else(subj_circSD_ES > 100, 100, subj_circSD_ES)) %>%
mutate(subj_circSD_ES = if_else(subj_circSD_ES < -100, -100, subj_circSD_ES)) %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_circSD_ES)) +
geom_histogram(binwidth = 5, fill = "red", alpha = 0.7) +
labs(x = "subj_circSD_ES [prior predictive distribution]",
subtitle = glue("p(< -100) = {sum(sim_datasets_unnest$subj_circSD_ES < -100)/length(sim_datasets_unnest$subj_circSD_ES)}\n p(>100) = {sum(sim_datasets_unnest$subj_circSD_ES > 100)/length(sim_datasets_unnest$subj_circSD_ES)}"))
plot_grid(subj_pre_circSD_plot, subj_post_circSD_plot, subj_circSD_ES_plot, ncol = 1)
narrower variance prior has reduced extreme values, good. might need to do more.
Calculate min and max subject pre, post, and circSD effect size in each simulation.
sim_subj_extremes <- sim_datasets %>%
unnest(dataset) %>%
mutate(subj_alpha1 = subj_alpha0 + subj_alphaD,
subj_alpha0_delinked = exp(subj_alpha0),
subj_alpha1_delinked = exp(subj_alpha1),
subj_circSD_effect = subj_alpha1_delinked - subj_alpha0_delinked ) %>%
group_by(sim_num) %>%
summarize_at(vars(subj_alpha0_delinked, subj_alpha1_delinked, subj_circSD_effect), list(max = max, min = min))
simSD <- sim_datasets %>%
unnest(dataset) %>%
mutate(subj_pre_circSD = exp(subj_alpha0),
subj_post_circSD = exp(subj_alpha0 + subj_alphaD)) %>%
group_by(sim_num) %>%
summarise_at(vars(subj_pre_circSD, subj_post_circSD), list(mean = mean, sd = sd))
plot_grid(
ggplot(simSD, aes(subj_pre_circSD_sd)) +
geom_histogram(binwidth = 2) +
xlim(0, 200) +
labs(x = "pre condition: observed SD of distribution of circSD, per sim",
subtitle = glue::glue("p(>50) = {mean(simSD$subj_pre_circSD_sd > 50)}")),
ggplot(simSD, aes(subj_post_circSD_sd)) +
geom_histogram(binwidth = 2) +
xlim(0, 200) +
labs(x = "post condition: observed SD of distribution of circSD, per sim",
subtitle = glue::glue("p(>50) = {mean(simSD$subj_post_circSD_sd > 50)}")),
ncol = 1
)
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 5 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
plot_grid(
ggplot(simSD, aes(subj_pre_circSD_mean, subj_pre_circSD_sd)) +
geom_point() +
xlim(0, 200) + ylim(0, 100) +
labs(x = "pre condition: observed mean of distribution of circSD, per sim",
y = "pre condition: observed sd"),
ggplot(simSD, aes(subj_post_circSD_mean, subj_post_circSD_sd)) +
geom_point() +
xlim(0, 200) + ylim(0, 100) +
labs(x = "post condition: observed mean of distribution of circSD, per sim",
y = "post condition: observed sd"),
ncol = 1
)
## Warning: Removed 11 rows containing missing values (geom_point).
## Warning: Removed 67 rows containing missing values (geom_point).
#What is the max subject pre condition value in each dataset
pre_plot <- sim_subj_extremes %>%
mutate(subj_alpha0_delinked_max = if_else(subj_alpha0_delinked_max > 120, 120, subj_alpha0_delinked_max)) %>%
ggplot() +
geom_histogram(aes(x = subj_alpha0_delinked_max, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "pre condition: max subj circSD per sim [subj_alpha0_delinked_max]",
subtitle = glue("p(<5) = {sum(sim_subj_extremes$subj_alpha0_delinked_max < 5)/length(sim_subj_extremes$subj_alpha0_delinked_max)}\n p(>120) = {sum(sim_subj_extremes$subj_alpha0_delinked_max > 120)/length(sim_subj_extremes$subj_alpha0_delinked_max)}"))
#What is the max subject post condition value in each dataset
post_plot <- sim_subj_extremes %>%
mutate(subj_alpha1_delinked_max = if_else(subj_alpha1_delinked_max > 120, 120, subj_alpha1_delinked_max)) %>%
ggplot() +
geom_histogram(aes(x = subj_alpha1_delinked_max, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "post condition: max subj circSD per sim [subj_alpha1_delinked_max]",
subtitle = glue("p(<5) = {sum(sim_subj_extremes$subj_alpha1_delinked_max < 5)/length(sim_subj_extremes$subj_alpha1_delinked_max)}\n p(>120) = {sum(sim_subj_extremes$subj_alpha1_delinked_max > 120)/length(sim_subj_extremes$subj_alpha1_delinked_max)}"))
plot_grid(pre_plot, post_plot, ncol =1, align = 'v')
narrower variance prior has reduced extreme values, good. might need to do more.
pre_plot <- sim_subj_extremes %>%
mutate(subj_alpha0_delinked_min = if_else(subj_alpha0_delinked_min > 120, 120, subj_alpha0_delinked_min)) %>%
ggplot() +
geom_histogram(aes(x = subj_alpha0_delinked_min, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "pre condition: min subj circSD per sim [subj_alpha0_delinked_min]",
subtitle = glue("p(<5) = {sum(sim_subj_extremes$subj_alpha0_delinked_min < 5)/length(sim_subj_extremes$subj_alpha0_delinked_min)}\n p(>120) = {sum(sim_subj_extremes$subj_alpha0_delinked_min > 120)/length(sim_subj_extremes$subj_alpha0_delinked_min)}"))
#What is the most extreme subject post condition value in each dataset
post_plot <- sim_subj_extremes %>%
mutate(subj_alpha1_delinked_min = if_else(subj_alpha1_delinked_min > 120, 120, subj_alpha1_delinked_min)) %>%
ggplot() +
geom_histogram(aes(x = subj_alpha1_delinked_min, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "post condition: min subj circSD per sim [subj_alpha1_delinked_min]",
subtitle = glue("p(<5) = {sum(sim_subj_extremes$subj_alpha1_delinked_min < 5)/length(sim_subj_extremes$subj_alpha1_delinked_min)}\n p(>120) = {sum(sim_subj_extremes$subj_alpha1_delinked_min > 120)/length(sim_subj_extremes$subj_alpha1_delinked_min)}"))
plot_grid(pre_plot, post_plot, ncol =1, align = 'v')
narrower variance prior has reduced extreme values, good. might need to do more.
#What is the most extreme change from pre to post in a subject
min_plot <- sim_subj_extremes %>%
mutate(subj_circSD_effect_min = if_else(subj_circSD_effect_min > 100 , 100, subj_circSD_effect_min)) %>%
mutate(subj_circSD_effect_min = if_else(subj_circSD_effect_min < -100 , -100, subj_circSD_effect_min)) %>%
ggplot() +
geom_histogram(aes(x = subj_circSD_effect_min, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "min subj circSD effect per sim \n[min(subj_alpha1_delinked - subj_alpha0_delinked)]",
subtitle = glue("p(<-100) = {sum(sim_subj_extremes$subj_circSD_effect_min < -100)/length(sim_subj_extremes$subj_circSD_effect_min)}\n p(>100) = {sum(sim_subj_extremes$subj_circSD_effect_min > 100)/length(sim_subj_extremes$subj_circSD_effect_min)}"))
#What is the most extreme subject post condition value in each dataset
max_plot <- sim_subj_extremes %>%
mutate(subj_circSD_effect_max = if_else(subj_circSD_effect_max > 100 , 100, subj_circSD_effect_max)) %>%
mutate(subj_circSD_effect_max = if_else(subj_circSD_effect_max < -100 , -100, subj_circSD_effect_max)) %>%
ggplot() +
geom_histogram(aes(x = subj_circSD_effect_max, fill = "red", alpha = 0.3), binwidth = 5) +
theme(legend.position = "none") +
labs(x = "max subj circSD effect per sim \n[max(subj_alpha1_delinked - subj_alpha0_delinked)]",
subtitle = glue("p(<-100) = {sum(sim_subj_extremes$subj_circSD_effect_max < -100)/length(sim_subj_extremes$subj_circSD_effect_max)}\n p(>100) = {sum(sim_subj_extremes$subj_circSD_effect_max > 100)/length(sim_subj_extremes$subj_circSD_effect_max)}"))
plot_grid(min_plot, max_plot, ncol =1, align = 'v')
narrower variance prior has reduced extreme values, good. might need to do more.
Both the max pre/post subject distributions of circSD and the max/min subject effect size distributions on circSD contained too many extreme values.
Can I see which prior distributions I need to change in order to correct these by plotting the points from those distributions along with their respective prior distribution samples?
(Is this a useful technique in general?)
Plots:
for (max pre circSD, max post circSD, max circSD ES, min circSD ES) * vs alpha0_mu * vs alpha0_sigma * vs alphaD_mu * vs alphaD_sigma
sim_datasets_join_extremes <- sim_datasets %>%
mutate(alpha0_mu_delinked = exp(alpha0_mu),
alpha1_mu_delinked = exp(alpha0_mu + alphaD_mu)) %>%
left_join(sim_subj_extremes, by = "sim_num") %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup()
#parameters describing uncertanity over group circ SDs in the pre condition (alpha0_mu, alpha0_sigma)
#vs
#implied pre condition min and max circSDs
preMax_vs_alpha0_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha0_delinked_max, y = alpha0_mu_delinked)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
preMin_vs_alpha0_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha0_delinked_min, y = alpha0_mu_delinked)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
preMax_vs_alpha0_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha0_delinked_max, y = alpha0_sigma)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
preMin_vs_alpha0_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha0_delinked_min, y = alpha0_sigma)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
plot_grid(preMin_vs_alpha0_mu, preMax_vs_alpha0_mu, preMin_vs_alpha0_sigma, preMax_vs_alpha0_sigma, nrow = 2, align = "hv")
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 69 rows containing missing values (geom_point).
## Warning: Removed 1 rows containing missing values (geom_point).
## Warning: Removed 69 rows containing missing values (geom_point).
#parameters describing uncertanity over change from pre to post condition (alphaD_mu, alphaD_sigma)
#vs
#implied post condition min and max circSDs
postMin_vs_alphaD_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha1_delinked_min, y = alphaD_mu)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
postMax_vs_alphaD_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha1_delinked_max, y = alphaD_mu)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
postMin_vs_alphaD_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha1_delinked_min, y = alphaD_sigma)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
postMax_vs_alphaD_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_alpha1_delinked_max, y = alphaD_sigma)) +
geom_point() +
scale_x_continuous(limits = c(0, 200))
plot_grid(postMin_vs_alphaD_mu, postMax_vs_alphaD_mu, postMin_vs_alphaD_sigma, postMax_vs_alphaD_sigma, nrow = 2, align = "hv")
## Warning: Removed 9 rows containing missing values (geom_point).
## Warning: Removed 221 rows containing missing values (geom_point).
## Warning: Removed 9 rows containing missing values (geom_point).
## Warning: Removed 221 rows containing missing values (geom_point).
#parameters describing uncertanity over group circ SDs in the pre condition (alpha0_mu, alpha0_sigma)
#vs
#implied effect size min and max
ESmin_vs_alpha0_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_min, y = alpha0_mu_delinked)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmax_vs_alpha0_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_max, y = alpha0_mu_delinked)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmin_vs_alpha0_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_min, y = alpha0_sigma)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmax_vs_alpha0_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_max, y = alpha0_sigma)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
plot_grid(ESmin_vs_alpha0_mu, ESmax_vs_alpha0_mu, ESmin_vs_alpha0_sigma, ESmax_vs_alpha0_sigma, nrow = 2, align = "hv")
## Warning: Removed 35 rows containing missing values (geom_point).
## Warning: Removed 221 rows containing missing values (geom_point).
## Warning: Removed 35 rows containing missing values (geom_point).
## Warning: Removed 221 rows containing missing values (geom_point).
#parameters describing uncertanity over change from pre to post condition(alphaD_mu, alphaD_sigma)
#vs
#implied effect size min and max
ESmin_vs_alphaD_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_min, y = alphaD_mu)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmax_vs_alphaD_mu <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_max, y = alphaD_mu)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmin_vs_alphaD_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_min, y = alphaD_sigma)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
ESmax_vs_alphaD_sigma <-
ggplot(sim_datasets_join_extremes, aes(x = subj_circSD_effect_max, y = alphaD_sigma)) +
geom_point() +
scale_x_continuous(limits = c(-100, 100))
plot_grid(ESmin_vs_alphaD_mu, ESmax_vs_alphaD_mu, ESmin_vs_alphaD_sigma, ESmax_vs_alphaD_sigma, nrow = 2, align = "hv")
## Warning: Removed 35 rows containing missing values (geom_point).
## Warning: Removed 221 rows containing missing values (geom_point).
## Warning: Removed 35 rows containing missing values (geom_point).
## Warning: Removed 221 rows containing missing values (geom_point).
Takeaways:
Unsurprisingly, effect size and pre,post extreme values would likely be better contained with narrower priors.
First ideas:
*Only narrow priors on mu parameters, alpha0_mu and alphaD_mu. This should help but the next step will likely be ncessary?
*Also narrow priors on the SD parameters, alpha0_sigma and alphaD_sigma.
(pre group mean, post group mean, post-pre mean change)
## linked prior distribution
plot_grid(
sim_datasets %>%
ggplot(aes(x = beta0_mu)) +
geom_histogram(binwidth = 0.05, fill = "#F16C66") +
theme(legend.position = "none",
axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
labs(x = "linked prior on group mean for pMem, pre condition",
title = "pMem group mean prior, beta0_mu") +
scale_x_continuous(limits = c(-5,5))
,
sim_datasets %>%
transmute(beta1_mu = beta0_mu + betaD_mu) %>%
ggplot(aes(x = beta1_mu)) +
geom_histogram(binwidth = 0.05, fill = "#685369") +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
labs(x = "linked prior on group mean for pMem, post condition",
title = "pMem group mean prior, beta0_mu + betaD_mu") +
scale_x_continuous(limits = c(-5,5))
,
nrow = 1,
align = "h"
)
## Warning: Removed 1 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 5 rows containing non-finite values (stat_bin).
## Warning: Removed 2 rows containing missing values (geom_bar).
##
## inv_logit(beta0_mu)
##
##
## inv_logit(beta0_mu + betaD_mu) = inv_logit(beta1_mu)
##
## delinked prior alpha0_mu
beta0_mu_delinked_stats <- sim_datasets %>%
transmute(delinked = exp(beta0_mu)/(exp(beta0_mu) + 1)) %>%
summarise(less_than_5 = sum(delinked < 0.05)/n(), greater_than_95 = sum(delinked > 0.95)/n())
plot_grid(
sim_datasets %>%
transmute(delinked = exp(beta0_mu)/(exp(beta0_mu) + 1)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 0.03, fill = "#F16C66") +
geom_vline(xintercept = c(0.05, 0.95), linetype = "dashed") +
labs(x = "prior on group mean for pMem, pre condition",
title = "pMem group mean prior, inv_logit(beta0_mu)",
subtitle = glue("prob(pMem < 0.05) = {beta0_mu_delinked_stats$less_than_5}\nprob(pMem > 0.95) = {beta0_mu_delinked_stats$greater_than_95}")) +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12))
,
sim_datasets %>%
transmute(delinked = exp(beta0_mu + betaD_mu)/(exp(beta0_mu + betaD_mu) + 1)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 0.03, fill = "#685369") +
labs(x = "prior on group mean for pMem, post condition",
title = "pMem group mean prior, inv_logit(alpha0_mu + alphaD_mu)") +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12))
,
nrow = 1,
align = "h"
)
##
## inv_logit(beta0_mu + betaD_mu) - inv_logit(beta0_mu) plot
##
betaD_mu_delinked_stats <- sim_datasets %>%
transmute(delinked = exp(beta0_mu + betaD_mu)/(exp(beta0_mu + betaD_mu) + 1) - exp(beta0_mu)/(exp(beta0_mu) + 1)) %>%
summarise(less_than_n80 = sum(delinked < -0.8)/n(), greater_than_80 = sum(delinked > 0.8)/n())
sim_datasets %>%
transmute(delinked = exp(beta0_mu + betaD_mu)/(exp(beta0_mu + betaD_mu) + 1) - exp(beta0_mu)/(exp(beta0_mu) + 1)) %>%
ggplot(aes(x = delinked)) +
geom_histogram(binwidth = 0.02, fill = "#00AEB2") +
geom_vline(xintercept = c(-0.8, 0.8), linetype = "dashed") +
labs(x = "prior on group mean for delta-pMem, mean post minus mean pre",
title = "delta-pMem group mean prior, inv_logit(beta0_mu + betaD_mu) - inv_logit(beta0_mu)",
subtitle = glue("prob(delta-pMem < -0.8) = {betaD_mu_delinked_stats$less_than_n80}\nprob(delta-pMem > 0.8) = {betaD_mu_delinked_stats$greater_than_80}")) +
theme(axis.title = element_text(size = 12),
plot.title = element_text(size = 12)) +
scale_x_continuous(breaks= pretty_breaks(10))
These priors still look good. no change from m1 to mean priors
(based on simulated subjects from each dataset)
(ignoring simulation sample sizes)
sim_datasets_unnest <- sim_datasets %>%
unnest(dataset) %>%
mutate(beta0_mu_delinked = exp(beta0_mu)/(exp(beta0_mu) + 1),
beta1_mu_delinked = exp(beta0_mu + betaD_mu)/(exp(beta0_mu + betaD_mu) + 1),
subj_pre_pMem = exp(subj_beta0)/(exp(subj_beta0) + 1),
subj_post_pMem = exp(subj_betaD + subj_beta0)/(exp(subj_betaD + subj_beta0) + 1),
subj_pMem_ES = subj_post_pMem - subj_pre_pMem)
pMem_ES_quantiles <- quantile(sim_datasets_unnest$subj_pMem_ES, c(0.025, 0.975))
plot_grid(
sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_pre_pMem)) +
geom_histogram(binwidth = 0.01, fill = "red", alpha = 0.3) +
labs(x = "subj_pre_pMem [prior predictive distribution]",
subtitle = glue("p(< 0.05) = {sum(sim_datasets_unnest$subj_pre_pMem < 0.05)/length(sim_datasets_unnest$subj_pre_pMem)}\n p(> 0.95) = {sum(sim_datasets_unnest$subj_pre_pMem > 0.95)/length(sim_datasets_unnest$subj_pre_pMem)}"))
,
sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_post_pMem)) +
geom_histogram(binwidth = 0.01, fill = "red", alpha = 0.3) +
labs(x = "subj_post_pMem [prior predictive distribution]",
subtitle = glue("p(< 0.05) = {sum(sim_datasets_unnest$subj_post_pMem < 0.05)/length(sim_datasets_unnest$subj_post_pMem)}\n p(> 0.95) = {sum(sim_datasets_unnest$subj_post_pMem > 0.95)/length(sim_datasets_unnest$subj_post_pMem)}"))
,
sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(x = subj_pMem_ES)) +
geom_histogram(binwidth = 0.01, fill = "red", alpha = 0.7) +
geom_vline(xintercept = pMem_ES_quantiles, linetype = "dashed") +
labs(x = "subj_pMem_ES [prior predictive distribution], (5%, 95%) quantile lines",
subtitle = glue("p(< -0.8) = {sum(sim_datasets_unnest$subj_pMem_ES < -0.8)/length(sim_datasets_unnest$subj_pMem_ES)}\n p(> 0.8) = {sum(sim_datasets_unnest$subj_pMem_ES > 0.80)/length(sim_datasets_unnest$subj_pMem_ES)}"))
,
ncol = 1
)
This also looks the same as m1.
Calculate min and max subject pre, post, and pMem effect size in each simulation.
sim_subj_extremes <- sim_datasets %>%
unnest(dataset) %>%
mutate(subj_beta1 = subj_beta0 + subj_betaD,
subj_beta0_delinked = exp(subj_beta0)/(exp(subj_beta0) + 1),
subj_beta1_delinked = exp(subj_beta0 + subj_betaD)/(exp(subj_beta0 + subj_betaD) + 1),
subj_pMem_effect = subj_beta1_delinked - subj_beta0_delinked ) %>%
group_by(sim_num) %>%
summarize_at(vars(subj_beta0_delinked, subj_beta1_delinked, subj_pMem_effect), list(max = max, min = min))
plot_grid(
#What is the max subject pre condition value in each dataset
sim_subj_extremes %>%
ggplot() +
geom_histogram(aes(x = subj_beta0_delinked_max, fill = "red", alpha = 0.3), binwidth = 0.02) +
theme(legend.position = "none") +
labs(x = "pre condition: max subj pMem per sim [subj_beta0_delinked_max]",
subtitle = glue("p(< 0.05) = {sum(sim_subj_extremes$subj_beta0_delinked_max < 0.05)/length(sim_subj_extremes$subj_beta0_delinked_max)}\n p(> 0.95) = {sum(sim_subj_extremes$subj_beta0_delinked_max > 0.95)/length(sim_subj_extremes$subj_beta0_delinked_max)}"))
,
#What is the max subject post condition value in each dataset
sim_subj_extremes %>%
ggplot() +
geom_histogram(aes(x = subj_beta1_delinked_max, fill = "red", alpha = 0.3), binwidth = 0.02) +
theme(legend.position = "none") +
labs(x = "post condition: max subj pMem per sim [subj_beta1_delinked_max]",
subtitle = glue("p(< 0.05) = {sum(sim_subj_extremes$subj_beta1_delinked_max < 0.05)/length(sim_subj_extremes$subj_beta1_delinked_max)}\n p(>0.95) = {sum(sim_subj_extremes$subj_beta1_delinked_max > 0.95)/length(sim_subj_extremes$subj_beta1_delinked_max)}"))
,
ncol =1,
align = 'v'
)
This slight skew towards 1 is better than the complete skew in m1 perhaps. This is a result of reducing variance.
plot_grid(
#What is the min subject pre condition value in each dataset
sim_subj_extremes %>%
ggplot() +
geom_histogram(aes(x = subj_beta0_delinked_min, fill = "red", alpha = 0.3), binwidth = 0.02) +
theme(legend.position = "none") +
labs(x = "pre condition: min subj pMem per sim [subj_beta0_delinked_min]",
subtitle = glue("p(< 0.05) = {sum(sim_subj_extremes$subj_beta0_delinked_min < 0.05)/length(sim_subj_extremes$subj_beta0_delinked_min)}\n p(> 0.95) = {sum(sim_subj_extremes$subj_beta0_delinked_min > 0.95)/length(sim_subj_extremes$subj_beta0_delinked_min)}"))
,
#What is the min subject post condition value in each dataset
sim_subj_extremes %>%
ggplot() +
geom_histogram(aes(x = subj_beta1_delinked_min, fill = "red", alpha = 0.3), binwidth = 0.02) +
theme(legend.position = "none") +
labs(x = "post condition: min subj pMem per sim [subj_beta1_delinked_min]",
subtitle = glue("p(< 0.05) = {sum(sim_subj_extremes$subj_beta1_delinked_min < 0.05)/length(sim_subj_extremes$subj_beta1_delinked_min)}\n p(>0.95) = {sum(sim_subj_extremes$subj_beta1_delinked_min > 0.95)/length(sim_subj_extremes$subj_beta1_delinked_min)}"))
,
ncol =1,
align = 'v'
)
This slight skew towards 0 is better than the complete skew in m1 perhaps. This is a result of reducing variance.
#What is the most extreme change from pre to post in a subject
min_plot <- sim_subj_extremes %>%
mutate(subj_pMem_effect_min = if_else(subj_pMem_effect_min > 1 , 1, subj_pMem_effect_min)) %>%
mutate(subj_pMem_effect_min = if_else(subj_pMem_effect_min < -1 , -1, subj_pMem_effect_min)) %>%
ggplot() +
geom_histogram(aes(x = subj_pMem_effect_min, fill = "red", alpha = 0.3), binwidth = 0.03) +
theme(legend.position = "none") +
labs(x = "min subj pMem effect per sim \n[min(subj_beta1_delinked - subj_beta0_delinked)]",
subtitle = glue("p(<-0.8) = {sum(sim_subj_extremes$subj_pMem_effect_min < -0.8)/length(sim_subj_extremes$subj_pMem_effect_min)}\n p(>0.8) = {sum(sim_subj_extremes$subj_pMem_effect_min > 0.8)/length(sim_subj_extremes$subj_pMem_effect_min)}")) +
xlim(c(-1, 1))
#What is the most extreme subject post condition value in each dataset
max_plot <- sim_subj_extremes %>%
mutate(subj_pMem_effect_max = if_else(subj_pMem_effect_max > 1 , 1, subj_pMem_effect_max)) %>%
mutate(subj_pMem_effect_max = if_else(subj_pMem_effect_max < -1 , -1, subj_pMem_effect_max)) %>%
ggplot() +
geom_histogram(aes(x = subj_pMem_effect_max, fill = "red", alpha = 0.3), binwidth = 0.03) +
theme(legend.position = "none") +
labs(x = "max subj pMem effect per sim \n[max(subj_beta1_delinked - subj_beta0_delinked)]",
subtitle = glue("p(<-0.8) = {sum(sim_subj_extremes$subj_pMem_effect_max < -0.8)/length(sim_subj_extremes$subj_pMem_effect_max)}\n p(>0.8) = {sum(sim_subj_extremes$subj_pMem_effect_max > 0.8)/length(sim_subj_extremes$subj_pMem_effect_max)}")) +
xlim(c(-1, 1))
plot_grid(min_plot, max_plot, ncol =1, align = 'v')
## Warning: Removed 2 rows containing missing values (geom_bar).
## Warning: Removed 2 rows containing missing values (geom_bar).
extreme effects are also contained within a more resonable range compared to m1
joint_pre_plot <- sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
ggplot(aes(y = exp(subj_alpha0), x = exp(subj_beta0)/(exp(subj_beta0) + 1))) +
geom_point() +
labs(x = 'pMem, pre', y = 'circSD, pre')
plot_grid(
joint_pre_plot,
ncol = 1,
align = 'v'
)
joint_ES_plot <- sim_datasets_unnest %>%
group_by(sim_num) %>%
sample_n(1) %>%
ungroup() %>%
mutate(subj_pMem_ES = exp(subj_beta0 + subj_betaD)/(exp(subj_beta0 + subj_betaD) + 1) - exp(subj_beta0)/(exp(subj_beta0) + 1),
subj_circSD_ES = exp(subj_alpha0 + subj_alphaD) - exp(subj_alpha0)) %>%
ggplot(aes(x = subj_pMem_ES, y = subj_circSD_ES)) +
geom_point() +
geom_vline(xintercept = 0, linetype = 'dashed') +
geom_hline(yintercept = 0, linetype = 'dashed')
plot_grid(
joint_ES_plot,
joint_ES_plot + ylim(c(-200, 200)) + labs(subtitle = "zoom in"),
ncol = 1,
align = 'v'
)
## Warning: Removed 14 rows containing missing values (geom_point).
(based on simulated obs from each dataset, ignoring subject labels)
sim_subj_obs_hist_count <- function(dataset, condition = 0){
dataset_obs <- dataset %>%
unnest(subj_obs) %>%
filter(stimulation == condition) %>%
select(obs_degree)
breaks <- seq(-180, 180, 5)
bincount <- hist(dataset_obs$obs_degree, breaks = breaks, plot = FALSE)$counts
bincount_names <- glue("c{breaks[-1]}")
names(bincount) <- bincount_names
bincount_df <- data.frame(as.list(bincount))
return(bincount_df)
}
make_quantmat <- function(sim_datasets, condition = 0){
bincounts <- sim_datasets %>%
select(dataset) %>%
mutate(subj_hist_counts = map(dataset, sim_subj_obs_hist_count, condition)) %>%
select(-dataset) %>%
unnest(subj_hist_counts) %>%
as_tibble()
xvals <- seq(-177.5, 177.5, 5)
probs <- seq(0.1,0.9,0.1)
quantmat <- as.data.frame(matrix(NA, nrow=ncol(bincounts), ncol=length(probs)))
names(quantmat) <- paste0("p",probs)
quantmat <- cbind(quantmat, xvals)
for (iQuant in 1:length(probs)){
quantmat[,paste0("p",probs[iQuant])] <- as.numeric(summarise_all(bincounts, ~quantile(., probs[iQuant])))
}
return(quantmat)
}
# calculate quantile mats from each condition
quantmat_cond0 <- make_quantmat(sim_datasets, 0)
quantmat_cond1 <- make_quantmat(sim_datasets, 1)
c_light <- "#DCBCBC"
c_light_highlight <- "#C79999"
c_mid <- "#B97C7C"
c_mid_highlight <- "#A25050"
c_dark <- "#8F2727"
c_dark_highlight <- "#7C0000"
plot_grid(
ggplot(quantmat_cond0, aes(x = xvals)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), color = c_dark, size = 1) +
scale_x_continuous(breaks=pretty_breaks(10)) +
labs(x = "error (degrees)", y = "count +/- quantile", subtitle = "without stimulation")
,
ggplot(quantmat_cond1, aes(x = xvals)) +
geom_ribbon(aes(ymax = p0.9, ymin = p0.1), fill = c_light) +
geom_ribbon(aes(ymax = p0.8, ymin = p0.2), fill = c_light_highlight) +
geom_ribbon(aes(ymax = p0.7, ymin = p0.3), fill = c_mid) +
geom_ribbon(aes(ymax = p0.6, ymin = p0.4), fill = c_mid_highlight) +
geom_line(aes(y = p0.5), color = c_dark, size = 1) +
scale_x_continuous(breaks=pretty_breaks(10)) +
labs(x = "error (degrees)", y = "count +/- quantile", subtitle = "with stimulation")
,
sim_datasets %>%
sample_n(1) %>%
unnest(dataset) %>%
unnest(subj_obs) %>%
filter(stimulation == sample(c(0,1), 1)) %T>%
print() %>%
ggplot(aes(x = obs_degree)) +
geom_density() +
labs(subtitle = "correctness check: random simulation, random condition") +
scale_x_continuous(breaks=pretty_breaks(10))
,
ncol = 1,
align = "v"
)
## # A tibble: 1,890 x 20
## sim_num alpha0_mu alpha0_sigma alphaD_mu alphaD_sigma beta0_mu
## <int> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 916 4.28 0.113 0.652 0.0992 -0.385
## 2 916 4.28 0.113 0.652 0.0992 -0.385
## 3 916 4.28 0.113 0.652 0.0992 -0.385
## 4 916 4.28 0.113 0.652 0.0992 -0.385
## 5 916 4.28 0.113 0.652 0.0992 -0.385
## 6 916 4.28 0.113 0.652 0.0992 -0.385
## 7 916 4.28 0.113 0.652 0.0992 -0.385
## 8 916 4.28 0.113 0.652 0.0992 -0.385
## 9 916 4.28 0.113 0.652 0.0992 -0.385
## 10 916 4.28 0.113 0.652 0.0992 -0.385
## # … with 1,880 more rows, and 14 more variables: beta0_sigma <dbl>,
## # betaD_mu <dbl>, betaD_sigma <dbl>, nsubj <dbl>, nobs_per_cond <dbl>,
## # subj <int>, subj_alpha0 <dbl>, subj_alphaD <dbl>, subj_beta0 <dbl>,
## # subj_betaD <dbl>, nobs_per_condition <dbl>, obs_radian <dbl>,
## # obs_degree <dbl>, stimulation <dbl>
looks good.